\documentclass{article}

\usepackage[french]{babel}
\usepackage[utf8]{inputenc}
\usepackage{textcomp}
\usepackage{amsmath}
\usepackage{float}
\usepackage[vlined,lined,linesnumbered,boxed,french]{algorithm2e}
\usepackage{graphicx}
\usepackage{stmaryrd}
\usepackage{upgreek}
\usepackage{listings}
\usepackage{subcaption}
\usepackage{subfig}
\usepackage[top=2cm, bottom=2.3cm]{geometry}

%%%%%%%%%%%%%%%% Lengths %%%%%%%%%%%%%%%%
\setlength{\textwidth}{15.5cm}
\setlength{\evensidemargin}{0.5cm}
\setlength{\oddsidemargin}{0.5cm}

%%%%%%%%%%%%%%%% Variables %%%%%%%%%%%%%%%%
\def\projet{5}
\def\titre{Méthodes d'interpolation et d'intégration / Interpolation de surfaces et splines cubiques}
\def\groupe{1}
\def\equipe{3}
\def\responsible{gucortes}
\def\secretary{gsimon}
\def\others{jfoucault, bsaintes, agautier}

\begin{document}

%%%%%%%%%%%%%%%% Header %%%%%%%%%%%%%%%%
\noindent\begin{minipage}{0.98\textwidth}
  \vskip 0mm
  \noindent
  { \begin{tabular}{p{7.5cm}}
      {\bfseries \sffamily
        Projet n°\projet} \\ 
      {\itshape \titre}
    \end{tabular}}
  \hfill 
  \fbox{\begin{tabular}{l}
      {~\hfill \bfseries \sffamily Groupe n°\groupe\ - Equipe n°\equipe
        \hfill~} \\[2mm] 
      Responsable : \responsible \\
      Secrétaire : \secretary \\
      Codeurs : \others
    \end{tabular}}
  \vskip 4mm ~

  ~~~\parbox{0.95\textwidth}{\small \textit{Résumé~:} \sffamily Le but de ce projet consiste en la réalisation d'un modèle représentant les flux d'air autour d'une aile d'avion, afin d'obtenir une carte des pressions autout de ladite aile.
    Dans un premier temps, nous traçerons la courbe de l'aile à partir de points donnés, par interpolation.
    Dans un second temps, nous développerons une méthode permettant de calculer la longueur d'une courbe. 
    Enfin, dans un dernier temps, nous utiliserons le profil laminaire de l'aile ainsi que sa longueur pour tracer la carte des pressions demandée.}
  \vskip 1mm ~
\end{minipage}

%%%%%%%%%%%%%%%% Main part %%%%%%%%%%%%%%%%
\section{Interpolation du modèle d'aile d'avion choisi}
Dans cette première partie nous avons implémenté une méthode d'interpolation par spline cubique. Cette méthode d'interpolation sera utilisée sur un profil d'aile \emph{NACA 747A315} créé par le \emph{National Advisory Committee for Aeronautics}, ancienne agence américaine chargée de recherches aéronautiques, depuis remplacée par la \emph{NASA}.
\subsection{Interpolation par splines cubiques}
  Les splines d’interpolations peuvent être de différents ordre, on les dira cubique lorsqu'ils sont d'ordre 3. Concrètement il s'agit de représenter la fonction caractérisant notre aile par un ensemble de polynômes de degrés 3 tels que :
  \begin{equation}
    \forall k\in [0,n], \forall t\in[x_k,x_{k+1}],\hspace{2em} f(t) = a\times t^3 + b\times t^2 + c\times t + d   \hspace{5em}   a,b,c,d\in\mathrm{I\!R}
  \end{equation}
  Ces splines sont ceux de plus petit degré qui permettent une approximation \emph{$C^2$} de plus elles ont de bonnes propriétés de régularité. C'est d'ailleurs en utilisant les conditions sur la continuité, la continuité de la dérivée première et celle de la dérivée seconde que l'on arrive à un système linéaire à $n+1$ inconnues et $n-1$ équations, où $n$ correspond au nombre de points données pour l'interpolation.\\ \\
  \indent En pratique, l'interpolation par spline cubique se fait en deux étapes: 
  \begin{itemize}
    \item On calcule les dérivées secondes en résolvant le système linéaire \ref{eq:syst-lineaire} 
    \item On détermine le polynôme grâce aux coefficients et aux dérivées précédemment déterminées.
  \end{itemize}
  Commençons par étudier le calcul des dérivées secondes.
  
  \subsubsection{Détermination des dérivées secondes}
    L'obtention des dérivées se fait par résolution du système suivant \ref{eq:syst-lineaire} 
     On a de plus : 
    \begin{itemize}
      \item $y_j$ : ordonnée du point j 
      \item $x_j$ : abscisse du point j
    \end{itemize}
    
\vspace{1em}
\indent Afin de pouvoir résoudre le système \ref{eq:syst-lineaire}, on pose $y"_0=y"_n=0$, ce qui nous donnera une spline cubique \emph{naturelle}.
Concrètement la résolution de ce système peut se faire de deux manières différentes.\\
    \indent La première consiste à reconstruire la matrice tridiagonale ainsi que les deux vecteurs colonnes puis à résoudre le système matriciel avec \emph{numpy.linalg.solve}. Cette fonction de \emph{Numpy} utilise une décomposition LU, se qui implique que cette fonction se fait en $\Theta(n^3)$ opérations.\\
    \indent La seconde solution consiste à utiliser l'algorithme fourni par le \emph{Numerical Recipes} dans le chapitre 3.3 et qui consiste à résoudre le système en \emph{"cascade"} de haut en bas puis de bas en haut en utilisant les deux conditions aux bornes que sont $y"_0$ et $y"_n$. La résolution du système se fait donc en temps linéaire en la taille du système, \emph{i.e}, en $\Theta(2\times n)$.
    \begin{equation}[H]
        \begin{pmatrix}
           \frac{x_2-x_0}{3} & \frac{x_2-x_1}{6} & 0 & \cdots & & \cdots & 0 \\
           \frac{x_1-x_0}{6} & \ddots & \ddots & \ddots & & & \vdots \\
           0 & \ddots & \ddots & \ddots & \ddots &  & \vdots \\
           \vdots & \ddots & \frac{x_j-x_{j-1}}{6} & \frac{x_{j+1}-x_{j-1}}{3} & \frac{x{j+1}-x_j}{6} & \ddots & \vdots \\
           \vdots & & \ddots & \ddots & \ddots & \ddots & 0 \\
           \vdots & & & \ddots & \ddots & \ddots & \frac{x_n-x_{n-1}}{6} \\
           0 & \cdots & & \cdots & 0 & \frac{x_{n-1}-x_{n-1}}{6} & \frac{x_n-x_{n-2}}{3} 
        \end{pmatrix}
        \begin{pmatrix}
           y"_1 \\
           \vdots \\
           y"_{j-1} \\
           y"_{j} \\
           y"_{j+1} \\
           \vdots \\
           y"_{n-1}
        \end{pmatrix}
        =
        \begin{pmatrix}
           \lambda_1 \\
           \vdots \\
            \\
           \lambda_j \\
           \\
           \\
           \lambda_{n-1}
        \end{pmatrix}
      \label{eq:syst-lineaire}
    \end{equation}
    \begin{equation}
      avec,\hspace{1em} \lambda_j = \frac{y_{j+1}-y_j}{x_{j+1}-x_j} - \frac{y_j-y_{j-1}}{x_j-x_{j-1}}
    \end{equation}
   

 \subsubsection{Création du polynôme interpolateur}
  Une fois les dérivées calculées, il est possible de déterminer pour chaque intervalle défini par deux points, de l'ensemble des points à interpoler, les coefficients du polynôme de degré 3 qui interpole au mieux notre la fonction, ici de l'aile, sur cet intervalle. Nous avons donc $n$ polynômes interpolateurs de la forme : 
  \begin{equation}
    f(x) = A\times x^3 + B\times x^2 + C\times x + D
  \end{equation}
  avec : 
  \begin{equation}
    A = \frac{y"_{j+1}-y"_j}{6.(x_{j+1}-x_j)}
  \end{equation}
  \begin{equation}
    B = \frac{x_{j+1}.y"_j - x_j.y"_{j+1}}{2x_{j+1}-2x_j}
  \end{equation}
  \begin{equation}
    C = \frac{-2.x_{j+1}^2.y"_j-x_{j+1}^2.y"_{j+1}-2.x_{j+1}.x_j.y"_j+2.x_{j+1}.x_j.y"_{j+1}+x_j^2.y"_j+2.x_j.y"_{j+1}-6.y_j+6.y_{j+1}}{6.(x_{j+1}-x_j)}
  \end{equation}
  \begin{equation}
    D = \frac{x_{j+1}^2.x_j.(2.y"_j+y"_{j+1})+a.(6.y_j-x_j^2.(y"_j+2.y"_{j+1}))-6.x_j.y_{j+1}}{6.(x_{j+1}-x_j)}
  \end{equation}

  La création des polynôme se fait en temps constant, donc la complexité de l'interpolation par spline cubique est $\Theta(n+2n)$, \emph{i.e}, $\Theta(n)$. \\
  \indent Pour ce qui est de la correction de notre fonction, la seule méthode mise en place ici, est graphique : on compare le tracé des points de la fonction à déterminer que l'on compare à celui de la fonction interpolé. La différence n'est pas flagrante dans la mesure ou l'espacement des points est suffisamment faible pour que leur seul tracé nous retourne quelque chose de suffisamment lisse. Nous avons donc également réalisé des tests, (voir figure \ref{fig:interpo}), en réduisant le nombre de points d'origine. Une autre méthode de test aurait été d'utiliser l'interpolation de \emph{Scipy} et de comparer les deux résultats.\\ \\

 \indent À l'origine nous utilisions des lambdas pour modéliser ces différents polynômes, l'objectif était de pouvoir utiliser ces lambda directement avec la fonction déterminant la longueur de courbe. Cependant à cause d'un problème de compréhension entre les membres du groupe, il s'est avéré que la dérivée de ces polynômes était nécessaire, or avec des \emph{lambda} obtenir la dérivée sans connaitre les coefficients est impossible. Ceci explique l'utilisation des lambda pour le tracé des différentes courbes, et l'utilisation des coefficients du polynôme pour le calcul de la longueur de courbe que nous étudierons plus tard dans ce rapport. L'ensemble de notre code aurait pu être uniformisé pour n'utiliser que les coefficients des polynômes.

  \begin{figure}[H]
    \begin{center}
    \includegraphics[width=0.48\textwidth]{../ressources/interpo.png}
    \includegraphics[width=0.48\textwidth]{../ressources/aile.png}
    \caption{Tracé des points de l'aile à droite, tracé de l'interpolation à gauche.}
   \label{fig:interpo}
    \end{center}
  \end{figure}

 

\section{Calcul de la longueur d'une courbe}

\subsection{Longueur de courbes}

\paragraph{}
Dans cette partie, nous abordons l'implémentation de méthodes mathématiques destinées à calculer la longueur de courbes représentant le profil d'une aile. Elles seront utilisées dans la partie suivante afin d'effectuer des calculs de divers paramètres physiques.

\paragraph{}
Une démonstration mathématique (qui sort du cadre de ce rapport) nous permet de montrer que la longueur d'une courbe décrite par la fonction $f$ entre les abscisses 0 et T se calcule de la façon suivante : 
\begin{equation}
L([0;T]) = \int ^T _0 \! \! \sqrt{1 + f'(x)^2}dx
\end{equation}
Les fonctions considérées étant polynomiales, trouver $f'$ est aisé; la difficulté réside donc dans le calcul de l'intégrale de la fonction $\sqrt{1 + f'(x)^2}dx$.

\subsection{Calcul numérique d'intégrales}

\paragraph{}
Les techniques de calcul d'intégrales sont nombreuses, mais toutes reposent sur le même principe : découper l'intervalle considéré en petits intervalles successifs de longueur appelée \textit{pas}, calculer une valeur approchée de l'intégrale sur chacun de ces intervalles, et sommer toutes ces valeurs. La partie importante consiste en le calcul de ces valeurs approchées, et différentes méthodes existent. \\
On peut citer comme exemple la plus connue : la méthode des rectangles gauche. Elle consiste à considérer la fonction $f$ constante sur les intervalles et ayant pour valeur la valeur de $f$ en l'éxtrémité gauche de l'intervalle. L'intégrale sur les intervalles est alors simple à calculer puisqu'il s'agit de l'aire du rectangle ayant pour cotés le pas et la valeur de la fonction.(cf figure \ref{fig:approx_rect_g})\\
\begin{figure}[!h]
\center
\includegraphics[width=3cm]{../ressources/approx_rect_g.png}
\caption{Exemple d'approximation de l'intégrale d'une fonction f sur un intervalle [a, b]}
\label{fig:approx_rect_g}
\end{figure}
\\
D'autres méthodes existent, telle que la méthode de Simpson, ayant pour avantage une meilleur approximation locale de la fonction au prix d'un coût plus élevé. Celui-ci n'est pas forcément négligeable car il est multiplié par le nombre de fois que la fonction est appelée, ce nombre augmentant lorsque la taille de l'intervalle augmente ou lorsque le pas diminue.

\paragraph{}
Nous avons étudié la convergence de ces méthodes en les testant sur différentes fonctions, et en faisant varier le pas. Les résultats obtenus pour la fonction $f(x) = \sqrt{x}$ sont résumés sur la figure \ref{fig:convergence_integrale_racine}. (d'autres courbes de convergence sont disponibles dans le dossier ressources)\\
\begin{figure}[!h]
\center
\includegraphics[width=8cm]{../ressources/racine.png}
\caption{Convergence des méthodes implémentées en fonction de la taille du pas ($f(x) = \sqrt{x}$)}
\label{fig:convergence_integrale_racine}
\end{figure}
\\
La première constatation importante est que toutes les méthodes considérées convergent vers la valeur exacte de l'intégrale. Cependant, cette convergence s'effectue à des vitesses différentes : en particulier, les méthodes de Boole-Villarceau et Simpson convergent rapidement, environ 20 subdivisions, là où les méthodes par rectangle demandent environ dix fois plus d'itérations.

\subsection{Découpage de l'intervalle}
Initialement le découpage de l'intervalle sur lequel on souhaitait intégrer la fonction était choisi arbitrairement. Nous avons perfectionné ce découpage suite à des idées présentées en cours. Désormais, une précision $\epsilon$ est fixée et on découpe l'intervalle en puissance de 2 croissantes tant que la différence en valeur absolue entre deux intégrales consécutives est supérieure à cet $\epsilon$. Cette méthode permet de limiter la découpe de l'intervalle à un minimum nécessaire au calcul, à une précision fixée, de l'intégrale; le tout sans connaitre la valeur exacte de celle-ci.


\subsection{Généricité}
Plutôt que d'imposer le choix de la méthode à l'utilisateur, nous avons choisi d'implémenter une fonction de calcul d'intégrale générique, prenant en paramètre une référence vers une fonction effectuant l'approximation. Ainsi, il est très simple de changer de méthode, et même d'en développer de nouvelles, car il suffit d'écrire une fonction effectuant le calcul désiré. Par exemple, la fonction utilisée dans le calcul par la méthode des rectangles gauche retourne $step * f(i)$, où $i$ désigne l'extremité gauche de l'intervalle considéré.

\paragraph{}
Maintenant que nous avons à notre disposition des moyens efficaces pour calculer des longueurs de courbes, nous pouvons progresser dans la modélisation du profil de l'aile.


\section{Caractéristiques du profil laminaire de l'aile d'avion}
\subsection{Écoulement laminaire autour de l'aile}
  Connaissant désormais une fonction interpolant l'aile, il est possible de tracer les lignes d'écoulement laminaire de l'air, en suivant la formule \ref{eq:ecoul-air}
  \begin{equation}
    g(x) = (1-\lambda)\times f(x) + \lambda\times h_max  \hspace{1em} \forall\lambda\in[0,1]
    \label{eq:ecoul-air}
  \end{equation}
  Enfin à partir des ces lignes d'écoulement de l'air, il est finalement possible en utilisant tous les outils définis précédemment de calculer la carte des pressions de l'aile étudiée.

\subsection{Cartes des pressions}
\indent La forme profilée de l'aile et l'écoulement de l'air sur cette aile implique que l'air proche de la surface de l'aile parcourra une plus grande distance que l'air plus éloigné. Or le temps pour parcourir cette distance sera le même quelque soit la position vertical par rapport à l'aile. Ainsi il est évident que l'air proche de l'aile aura une vitesse plus élevé que l'air plus éloigné. Enfin d'après la formule de la pression \ref{eq:pression}, la pression sera plus forte sur le dessus de l'aile. C'est de la différence entre la pression au-dessus de l'aile et celle du dessous que nait la force de portance.
\begin{equation}
  P = \frac{1}{2}\times\rho\times V^2
  \label{eq:pression}
\end{equation}

Grâce aux outils implementés précédemment le calcul de la pression pour chaque ligne d'air est possible. Il suffit alors de remplir une matrice contenant les informations de pression puis de l'afficher.\\
\indent Nous avons rencontré une importante difficulté au moment de joindre les différentes parties, notamment au niveau du format des données. Tout d'abord, les équations des courbes d'air sont fragmentées à cause du procédé d'interpolation. Il faut donc utiliser le bon polynôme d'interpolation sur le bon intervalle. Ensuite, les courbes d'écoulement de l'air sont contenues dans un plan de dimensions 1x1. Or, la matrice est composée d'un nombre entier de lignes et colonnes, il a fallu là encore accorder les indices de la matrice avec les cordonnées du plan de l'aile.\\
Nous n'avons donc pas anticipé les problèmes rencontrés et il en résulte un code peu propre, écrit à la limite de l'empirisme mais qui donne néanmoins des résultats très convenables. En effet, bien que le centre de l'aile aurait pu être remplie en noir afin de voir plus précisément les forces de pression s'exerçant sur l'aile, on remarque tout de même une différence de pression entre l'intrados et l'extrados de l'aile qui induit la force de portance comme le montre la figure \ref{fig:cart_pression}.
\begin{figure}[!h]
\center
\includegraphics[width=0.5\textwidth]{../ressources/carte_pressions.png}
\caption{Carte des pression de l'aile.}
\label{fig:cart_pression}
\end{figure}


\end{document}
